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We consider localisation problems belonging to the chiral symmetry classes, in which sublattice 
symmetry is responsible for singular behaviour at a band centre. We formulate models which have 
the relevant symmetries and which are generalisations of the network model introduced previously 
in the context of the integer quantum Hall plateau transition. We show that the generalisations 
required can be re-expressed as corresponding to the introduction of absorption and amplification 
into either the original network model, or the variants of it that represent disordered superconduc- 
tors. In addition, we demonstrate that by imposing appropriate constraints on disorder, a lattice 
version of the Dirac equation with a random vector potential can be obtained, as well as new types 
of critical behaviour. These models represent a convenient starting point for analytic discussions 
and computational studies, and we investigate in detail a two-dimensional example without time- 
reversal invariance. It exhibits both localised and critical phases, and band-centre singularities in 
the critical phase approach more closely in small systems the expected asymptotic form than in 
other known realisations of the symmetry class. 
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I. INTRODUCTION 

The classification by symmetry of Hamiltonians for dis- 
ordered systems provides an important framework in the 
study of Anderson localisation. Three standard sym- 
metry classes are long-established, and represented in 
the zero-dimensional limit|-|by the three Wigner-Dyson 
random matrix ensemblesB The existence of additional 
symmetry classes has been implicit in work which also 
has an extensive history, but a complete classification 
has been sttt out only rather recently, by Altland and 
Zirnbauerim Hamiltonians belonging to one of these ad- 
ditional symmetry classes are characterised by a discrete 
symmetry which relates eigenvalues and eigenfunctions 
in pairs. This paper is concerned with three examples, 
known as the chiral symmetry classes. They can be re- 
alised as tight-binding Hamiltonians for systems in which 
sites can be divided into two sublattices and non-zero 
matrix elements connect only sites from opposite sublat- 
tices. In these terms, the discrete symmetry operation is 
multiplication of wavefunction components on one sub- 
lattice by a factor of minus one: applied to an eigenstate 
with energy E, this generates a second eigenstate, with 
energy —E. Three distinct such symmetry classes are 
possible, termed!! AIII for systems without time-reversal 
invariance, BDI for spinless systems with time-reversal 
invariance, and CII for systems having time-reversal but 
not spin-rotation symmetry. 

Early work on disordered Hamiltonians with chiral 
symmetry followed a variety of independent directions. 
A one-dirnensional model of this type was solved exactly 
by Dysonta, who found a divergence in the density of 
states at the band centre, of a form shown by Ovchin- 
nikov and ErikhmanEl to be characteristic for the sym- 
metry class in one dimension. A closely related feature 
of this and similar models is that, while states away from 
the band centre are localised, both the mean free path 



and the localisation length diverge on approaching the 
band centre, a^revealed in early, treatments by Eggarter 
and RiedingerQand by Ziman.0 Two-dimensional sys- 
tems with- chiral symmetry were studied by Gade and 
Wegner,Btl using a mapping to a non- linear sigma model: 
their calculations for the class AIII with weak disorder 
indicate a strongly divergent density of states at the band 
centre, also accompanied by a divergence of the localisa- 
tion length. Chiral random—matrix ensembles, investi- 
gated tw Nagao and SlevinEHl and by Verbaarschot and 
Zahedjij likewise show band centre anomalies in the den- 
sity of states, albeit on energy scales comparable with 
the level spacing. And in a separate development, ex- 
tensive analytical results have been obtained by Ners- 
esyan et at3 and by Ludwig et aL3 for a particular two- 
dimensional model with chiral symmetry, the massless 
two-dimensional Dirac equation with random vector po- 
tential. 

Much subsequent work has helped to establish the 
range of behaviour possible in these systems. In one 
dimensional models, a finite localisation length at zero 
energy and a modification of the Dyson singularity in 
the density of states can be inducedElESEJ by terms in 
the Hamiltonian that in the absence of disorder would 
generate a spectral gap. For two-dimensional systems, it 
is recognised that a random vector potential in the Dirac 
equation constitutes a special choice of disorder, and scal- 
ing flow from this towards a generic fixed point, with 
properties similar to those derived by, Gade, has been 
studied in a field-theoretic frameworkli3 and by using re- 
sults on random Gaus sian ., su rfaces £3 Moreover, while 
early numerical resultsE3ii3u3 were quite different from 
analytical predictions, careful tuning of model parame- 
ters and study of large systems has recently brought cal , 
culations closer to expected asymptotic behaviour. tZIEil 
Numerical work has also shown that a finite value for the 
zero-energy localisation length can be produced in two- 
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dimensional chiral models, bji^the same mechanism that 
is effective in one dimension.!!!! 

Our aim in the following is to formulate and study 
new representatives of the pchiral symmetry classes, in 
the form of network modelsr3 In place of a Hamiltonian, 
these use the ideas of scattering theory and may be spec- 
ified by a transfer matrixJ23 or by a uniiacy. evolution op- 
erator for one step of discretised time.B3EJ The versions 
we set out here are connected in two distinct ways to par- 
ticular network models without chiral symmetry which 
have been studied previously. For the class AIII, these 
connections are to the U(l) network mode£3 investigated 
in the context of the integer quantum Hall plateau transi- 
tion, while for CII and BDI, respectively, they are to the 
SU(2) and 0(1) network models, which describe plateau 
transitions in dirty superconductors E3c3'E3 In each case, 
the model with chiral symmetry is constructed by cou- 
pling two copies of the partner model. This coupling 
can be re-expressed after a transformation as equivalent 
to introduction of absorption and coherent amplification 
in the original -models. This equivalence parallels the 
established linlcj between chiral symmetry classes and 
non-Hermitian random operators. 

The remainder of this paper is organised as follows. 
We show in Sec.|| how to construct network models be- 
longing to the chiral symmetry classes, and discuss in 
Sec. Ill some basic aspects, including the relation to non- 
Hermitian models, the continuum limit, and a lattice ver- 
sion of disorder analogous to a random vector potential 
in the Dirac equation. We present results from a numer- 
ical study of a two-dimensional model in the symmetry 
class AIII in Sec.|v|, and summarise in Sec.fv]. 



quirement that T. X HY, X = —H. Taking the scattering 
matrix to have the symmetry of e lH , this implies that 
T, X SY, X = S^ 1 . For the transfer matrix, it appears at 
first that one has a choice, taking flux in the pairs of 
scattering channels on which Y, x acts to propagate either 
in opposite directions or in the same direction. In fact 
only the former is tenable. It implies that the anticom- 
mutator {J,T, X } = and that the chiral symmetry con- 
dition is S^TEa; = T, a property preserved under matrix 
multiplication. By contrast, the latter choice would lead 
to the commutator [S a , T] = and the chiral symmetry 
condition T* X TY, X = T . Since the last condition is not 
in general preserved under matrix multiplication, connec- 
tion in series of two systems of this kind would generate 
a sample without chiral symmetry, and we therefore re- 
ject this alternative. Systems in class AIII have no other 
relevant discrete symmetry; those in classes BDI and CII 
are also invariant under time-reversal, in the absence and 
presence of Kramers degeneracy respectively. 



A. Network model for class AIII 

We recall first the essential features of the U(l) net- 
work model for the integer quantum Hall plateau transi- 
tion. A wavefunction in this model takes complex values 
zi on the links I of the lattice illustrated by the full lines 
of Fig. ^. The forms of the transfer matrix and of the evo- 
lution operator follow from the properties of the elemen- 
tary building units, shown in Fig. [lj A particle acquires 
a phase <f>i on traversing link /, so that in a stationary 
state amplitudes at either end are related by 



z' = z . 



II. CONSTRUCTION OF NETWORK MODELS 
WITH CHIRAL SYMMETRY 

In this section we formulate network models for each of 
the chiral symmetry classes, giving a detailed treatment 
of the chiral unitary class (AIII) , and indicating in outline 
the equivalent steps for the chiral orthogonal (BDI) and 
chiral symplectic (CII) classes. Our strategy is simply 
to construct the two-dimensional internal space associ- 
ated with chiral symmetry using two related copies of a 
network model without that symmetry. For each sym- 
metry class we focus on two-dimensional models; models 
in quasi-one and three dimensions can be constructed in 
the same way, starting for example from Eq. (0). 

The symmetry of a disordered system may be discussed 
in terms of a Hamiltonian H, a scattering matrix S, or 
a transfer matrix T. We are concerned with systems 
which conserve probability density, so that the Hamil- 
tonian is Hermitian and the scattering matrix is uni- 
tary: W = H and = S^ 1 . The equivalent condi- 
tion for the transfer matrix involves the current operator 
J and reads T _1 = JT^J. Chiral symmetry is imple- 
mented on a two-dimensional internal space in which the 
Pauli operator T, x acts: for the Hamiltonian it is the re- 



(1) 



In a similar way, stationary state amplitudes on the four 
links which meet at a node are related by a 2 x 2 transfer 
matrix 



cosh(a) sinh(a) 
sinh(a) cosh(a) 



(2) 



where a is real and all phase factors are associated with 
links. This equation may re-written in terms of a scat- 
tering matrix as 



cos(a) — sin(a) 
sin(a) cos(a) 



(3) 



with sin(a) = tanh(a). The transfer matrix that re- 
sults from assembling these units is described in de- 
tail in Ref. 22, and the time evolution operator in 
Refs.|| andgj. 

Introduction of a two-dimensional internal space asso- 
ciated with chiral symmetry results in a doubling of the 
number of wavefunction components. In this way, start- 
ing from two copies of the U(l) model the link ampli- 
tudes become two-component complex numbers, z/. In 
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FIG. 1: Units from which the network model is constructed: 
a link (left) and a node (right). 



place of Eq. ([!]), the scattering properties of a link are 
characterised by a 2 x 2 transfer matrix T, with 



z' = Tz . 

Requiring Y, x TYj x = T, T has the form 

T _ ( cosh(6) sinh(&) 
sinh(6) cosh(fe) 



(4) 



(5) 



where is a real phase and b is a real hyperbolic angle. 
It remains to discuss scattering at nodes of the doubled 
system. We replace Eq. (0) by 



zi 

Z2 



= 1 



cosh(a) sinh(a) 
sinh(a) cosh(a) 



Z.3 
Z4 



(6) 



where 1 is here the unit matrix in the two-component 
space introduced on links. This choice amounts to the 
most general one compatible with chiral symmetry, since 
all scattering within the two-component space may be 
included in the link transfer matrices, T. 

Combining these elements to make a two-dimensional 
system, we arrive at the model shown schematically in 
Fig. ^|. The transfer matrix for the system as a whole 
acts in the [1, 1] (or [1,T]) direction and may be written 
as a product of factors relating amplitudes on successive 
slices of the system. Alternate factors in the product 
represent links and nodes, and consist respectively of re- 
peated versions of the 2x2 and 4x4 blocks appearing 
in Eqs. (||) and (ft). We introduce disorder by taking the 
phase 4> in Eq. (a) to be an independent, uniformly dis- 
tributed random variable on each link, and take the pa- 
rameter a, characterising scattering at nodes, to be non- 
random. We have considered two ways to set the value 
of the coupling between the chiral subspaces, which we 
parameterise in terms of either the hyperbolic angle ap- 
pearing in Eq. (|J) or the compact angle (3 related to b 
by sin(/3) = — tanh(fr): we take (3 either uniformly dis- 
tributed, or b to have a normal distribution of variance 
g!E\ We ensure that the system is statistically invariant 
under ir/2 rotations of the lattice, which places require- 
ments on the node parameter o, exactly as in the U(l) 
model: nodes lie on two distinct sub-lattices, and the 



node parameter a on one sublattice is related to the pa- 
rameter a' on the other sublattice by the duality relation 
sinh(a) sinh(a') = 1. 



FIG. 2: Structure of the AIII network model. Full and dashed 
lines indicate the two U(l) models from which the system is 
constructed, with nodes located at the vertices of the two lat- 
tices. Scattering that couples the sub-systems is represented 
schematically by square boxes, shown only in the upper left 
plaquette. In one step of time evolution, flux propagates be- 
tween successive points marked with filled circles, in the di- 
rections indicated with arrows. 



The same system can alternatively be described using a 
time evolution operator, S, in place of a transfer matrix. 
In order to specify this operator, which is unitary and 
has the symmetry of a scattering matrix, it is convenient 
hrst to consider the special case b — 0, in which the 
two copies of the network model are uncoupled. Let U 
denote the time evolution operator for one copy. Then, 
since from Eq. (fj) link phases are the same in both copies 
but propagation directions are opposite, the evolution 
operator for the other copy is W . The dimension of U is 
equal to the number of links in the system. It is useful 
to define a diagonal matrix of the same dimension, with 
the angles Pi for each link I as diagonal entries: we use 
(3 to denote this matrix. The time-evolution operator for 
the system with chiral symmetry can then be written 



Ucos(p) 
sin(/3) 



-Z7sin(/3) W 
cos(/3) W 



It is straightforward to check that S^^E^ 
that Sft = S- 1 . 



(7) 



S' 1 and 
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B. Network model for classes BDI and CII 



A. Spectral symmetries 



A model in class BDI can be obtained from one in class 
All! simply by imposing time reversal invariance as an 
additional symmetry. This condition is conventionally 
written in the form H* = H, but for a discussion based 
on scattering matrices it is more convenient to make the 
transformation H — ► QHQ^ 1 with Q = exp(i?£a;). 
This transformation leaves the chiral symmetry relation 
T, x HYi x — —H unchanged. In the transformed basis one 
has H* = -H, S = S* and T* = T. To ensure a real 
time evolution operator, we restrict the link phases <f> to 
the values and tt. Choosing these values randomly, the 
BDI model consists of two coupled copies of the class 
D models studied-jrecently in the context of disordered 
superconductors.EZl Alternatively, one could set <j> = on 
all links, and introduce disorder only through the chiral 
couplings (3. 

We turn next to the symmetry class CII. Kramers de- 
generacy is a defining feature of the class and necessitates 
the introduction of an additional two-dimensional space 
arising from spin. The time reversal operation includes 
reversal of spin direction. Defining C = ir y , where r y is 
a Pauli matrix acting in the additional space, it is con- 
ventionally written in the form CH*C~ 1 =H. As for the 
class BDI, it is again convenient to make the transforma- 
tion H — > QHQ^ 1 with Q = exp (ijTi x ). In the trans- 
formed basis one has CH*^ 1 = -H, CSX- 1 = S and 
CT*C _1 = T as equivalent expressions of time-reversal 
invariance. Applying these ideas to a network model, 
four channels propagate on a single link which, generalis- 
ing Eq. (^) , has a 4 x 4 transfer matrix T with the generic 
form 



T = v. 



cosh(6) sinh(fo) 
sinh(6) cosh(fr) 



(8) 



where v is an SU(2) matrix and & is a real hyperbolic an- 
gle. Adopting this form, the time evolution operator for 
class CII has the structure given in Eq. tfjh, but with U 
representing a class C network model, studied previously 
in connection with the spin quantum Hall effect .E3c3 The 
links of a such a class C model carry two co-propagating 
channels, coming from two spin components, and the evo- 
lution operator satisfies CC/*C _1 = U. 



III. DISCUSSION OF THE MODELS 



The eigenphases E a and eigenvectors ip a of the time 
evolution operator S play the same role for a network 
model as do energy eigenvalues and eigenvectors for a 
model defined using a Hamiltonian. They satisfy 



Si/> a =e tB -il> a . 



(9) 



Chiral symmetry, Y, X ST, X = S^ 1 , has the consequence 
that T, x ip a is an eigenvector of S with opposite eigen- 
phase: 



ST, x ip 



(10) 



Because of this, the spectrum of eigenphases is symmet- 
ric around the point E = and their density p(E) may 
develop a singularity there. In addition, for models from 
the class CII, each eigenphase is Kramers degenerate. 
Specifically, in addition to ip a , S has a second eigenvec- 
tor, CE x ipa, with the same eigenvalue e lE ". 

The symmetries described in the preceding paragraph 
are the ones of fundamental interest in this work. How- 
ever, all network models based on the lattice of Fig.^ 
have in addition a bipartite structure. This structure 
leads to symmetry of the eigenphase spectrum under the 
translation E — > tt + E. When combined with chiral sym- 
metry, E — > —E, the consequence is that behaviour near 
E = is mirrored exactly at E — tt, and that similar 
behaviour is repeated at E = ±7r/2. |— . 

In detail, the bipartite structure (which is inherited^ 
from the network model for the quantum Hall plateau 
transition, and discussed further in Sec. HID) has the 
consequence that the evolution operator U appearing in 
Eq. (|7]) (or its counterparts for classes BDI and CII) can 
be written with a suitable ordering of the basis states as 



U 



A 
B 



(11) 



Introducing in the same basis a diagonal matrix a = 
diag(+l, — 1), bipartite lattice structure implies that 
alia = —U. Applying this to our models (and suitably 
extending the meaning of a), one finds that the pairs of 
eigenvectors of S identified above, ip a and 'E x ^> a , are re- 
lated to pairs aip a and aT, x ip a , with eigenvalues tt + E a 
and 7r — E a respectively. These relationships are illus- 
trated in Fig.^. They result in an eigenphase spectrum 
that is reflection-symmetric not only about the points 
E = and E = tt but also about the points E = ±tt/2. 



In this section we discuss some basic aspects which 
are common to models from all three chiral symmetry 
classes. These include: symmetry of the spectrum of the 
time-evolution operator; the relation to network mod- 
els which have absorption and amplification; the Green 
function; the continuum limit; behaviour for some spe- 
cial parameter values; and a lattice analogue of random 
vector potential disorder in the Dirac equation. 



B. Chiral symmetry and network models with 
absorption and amplification 

Any Hamiltonian H which satisifes the chiral sym- 
metry condition Y> x HYi x = —H has a block structure 
that is most clearly displayed after making the rotation 
R = {T, x + E^)/V2. In the rotated basis, the chiral 
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FIG. 3: Symmetries of the spectrum of the time evolution 
operator for a network model defined on the lattice illustrated 
in Pig Four eigenphases in the spectrum of S are shown 
as points on the unit circle in the complex plane. They are 
related by chiral symmetry combined with the symmetry of 
the bipartite lattice. 



symmetry condition involves S z rather than E x , since 
S z = KL X R- X . We have 



RHR- 1 



h 

ht 



(12) 



This establishes a decomposition of the Hermitian oper- 
ator H into the operators h and , which are in general 
non-Hermitian. We show in this section that for a net- 
work model with chiral symmetry there exists a similar 
decomposition into a pair of systems, each with half the 
number of degrees of freedom. A single member of the 
pair, taken in isolation, can be interpreted as a system in 
which there is absorption and amplification of flux. We 
discuss this separation first in the context of the transfer 
matrix, and then in terms of the time evolution operator. 

The logic can also be applied in the opposite directicm- 
In this case, just as Eq. ( |T2| ) provides a useful relationE2l 
between a non-Hermitian operator h, and the Hermitian 
operator H , so the arguments below can be used to relate 
a network model without flux conservation to a doubled 
system in which flux is conserved. 



1. Transfer matrices with absorption and amplification 

For our models all factors entering the transfer matrix 
are diagonal in the two-dimensional space on which Y, x 
acts, except for the 2x2 matrices appearing in Eqs. (|^) 
and (||). These matrices are diagonalised by the rota- 
tion R, which is disorder independent and which mixes 
counter-propagating channels: 



R 



cosh(6) sinh(6) 
sinh(6) cosh(6) 



R- 







(13) 



whole. The result reads symbolically 



RTR- 









T(-b) 



(14) 



In the case of our model with AIII symmetry, T(b) is 
obtained from the transfer matrix for the correspond- 
ing U(l) model by including imaginary parts in the link 
phases: 



ibi 



(15) 



In a similar way, for our models with BDI or CII symme- 
try, T(b) is obtained from the transfer matrices for class 
D or class C models respectively, by associating an extra 
factor of exp(6/) with propagation along each link I. 

If T(b) or T(—b) is regarded as the transfer matrix 
of a system in its own right, then that system is one 
with absorption and amplification. In fact, of course, 
the rotation R gives the current operator J off-diagonal 
components which couple T(b) and T(—b), so that T, 
properly considered, is current-conserving as it should 
be. 

From Eq. ( |l4| ) , one sees that the spectrum of Lyapunov 
exponents for a network model with chiral symmetry 
consists of separate contributions, from T(b) and from 
T{—b). Models in which the 6/ are random and dis- 
tributed symmetrically about bi = clearly have doubly 
degenerate Lyapunov exponents. These can be calcu- 
lated by studying just one member of the pair T(±6), 
which is important because it halves the size of matrices 
involved in numerical calculations. 



2. Time evolution operators with absorption and 
amplification 

It is interesting to re-examine this separation of the 
chiral model into subsystems, starting from the time- 
evolution operator in place of the transfer matrix. We 
begin by writing the action of the time-evolution opera- 
tor S in a way that emphasises the 2x2 block structure 
evident in Eq. (0) : 



Z, 

Z-2 



= s 



Zi 
Zi 



(16) 



The dimensions of the vectors Z and Z appearing in 
Eq. (|l^) are the same as those of the matrix U in Eq. (0). 
This expression can be reorganised into the form 



where K is 
K 



Z 3 
Zi 



U 
u 



= K 



Zi 
Z 2 



cosh(£>) sinh(6) 
sinh(6) cosh(6) 



(17) 



(18) 



Because this rotation is independent of disorder, we can 
apply it to the transfer matrix T for the system as a 



Properly speaking, since Z\ and Zi are amplitudes of 
oppositely propagating fluxes, K should be interpreted as 
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a transfer matrix, in the way illustrated in Fig. |]. Despite 
this, our approach is to apply to K the rotation R, and 
to discuss the rotated matrix as if it were an evolution 
operator. Doing this, we find 



where 



JlhI ' - { U(-b) 



U{±b) = Uexp{±b). 



(19) 



(20) 



For the special case 6 = 0, the subsystems represented 
by full and dashed lines in Fig. [2] are uncoupled. In this 
case one has U(Q) = U, which is the time evolution op- 
erator for one subsystem, and its inverse for the other 
subsystem. For general b, one can regard U(±b) as time 
evolution operators for systems with absorption and am- 
plification. However, it is only the stationary states of 
these systems which have a simple physical significance. 
Such states exist either for an open system with scat- 
tering boundary conditions (the situation represented by 
the transfer matrix T), or for a closed system if a pair 
of eigenphases are zero (which can be arranged by fine- 
tuning system parameters). In this case, we can define 
eigenvectors 



Z+=U(b)Z+ and Z_ =U{-b)Z_ 



(21) 



Undoing the rotation R, the stationary states of S satisfy 
tp± = Sip±, with 



and ip- 



-Z_ 



(22) 



u 



z, 



u 



FIG. 4: Schematic representation of the action of the time 
evolution operator S and the matrix K, appearing in Eqs. (^) 
and Scattering amplitudes within the chiral space are 

denoted by c = cos(/3) and s = sin(/3). 



C. Green function 

The spectral symmetry of the time evolution operator 
for network models with chiral symmetry is also shown 
clearly by considering the Green function, defined for 
complex energy z by 



G[z) 



1 

1 - zS 



(23) 



The combined consequences of the unitarity and chiral 
symmetry of S together fix the form of G(z), which is 
particularly simple at z = 1. For a discussion of thiSj it 
is again convenient to make the rotation R. We set S = 
RSR- 1 and G(z) = RG(z)R-\ so that T, Z ST, Z = ST 1 . 
Unitarity of S implies for G{z) 

G\z*)+G{z~ l ) = 1, (24) 

while chiral symmetry is the requirement that 

Y, z G{z)Y. z + G{z- 1 ) = l. (25) 

At z = 1, these two equations reduce to &(1) + G(l) = 1 
and S 2 G(1)E Z + G(l) = 1. Their solution fixes the form 
for G(l) to be 



(26) 



where D is a complex matrix. 

Our aim, then, is to determine G(l) explicitly for the 
models we have introduced. We start from Eq. fl22|), 
which we write as 



G{z) = zSG{z) + 1 . 

Following the step from Eq. (|l^) to Eq. (|l 
matrix K in place of S, and defining 



K 



z 
1 



we obtain 



G(z) 



1 



1 -K(z)' 



(27) 
using the 

(28) 
(29) 



This result simplifies further at z = 1 and after rotation 
byi? 



G(l) 



Note that, since 



1 

Ue- b 



1 - Ue- 



1 + Ue b 
1 - Ue b 
1 



(30) 



"1 + Ue b ~ 


t 


'l + Ue- b ~ 


1 - Ue b _ 




1 - Ue~ b _ 



(31) 



it is clear the Green function has the form expected from 
Eq. (H). 

Finally, we remark that similar manipulations may be 
used to construct a Hermitian operator which has the 
same eigenvectors as S and has the block structure 



U 



11-5 
2^1 + 5 



1 

27 





1 - Ue b 
1 + Ue b 



1 - Ue 



Ue~ b 




(32) 



7 



D. Continuum limit 

The network model representing the quantum Hall 
plateau transition is knowro to have a continuum limit 
in which the time evolution operator is that for a Dirac 
Hamiltonian in two space dimensions, with random vec- 
tor and scalar potentials and a uniform mass which varies 
through zero on crossing the transition. In this subsec- 
tion we investigate the continuum limit for network mod- 
els with chiral symmetry. As one might expect, this en- 
ables us to make links with previous work on thfi D irac 
equation with randomness in the classes AIIlEJ'EJ and 
BDI.E3 A similar approach should also be of interest for 
the CII network model, but we do not explore it here. 

We consider first the model for the class AIII. We ob- 
tain the continuum limit by taking all link phases and 
all scattering amplitudes within the chiral space to be 
small,y so that in Eq. (||) \4>\ < 1 and |b| < 1, and by 
taking the node parameter close to its self-dual value, so 
that in Eq. a = 7r/4 + m/2, with \m\ < 1. 

Before going further, it is necessary to deal with a tech- 
nical point which arises because of the bipart ite structure 
of the models, discussed above in Sec. [II A. This struc- 



ture has its origin in the form of U displayed in Eq. (11), 
but is obscured in S because of the nature of the off- 
diagonal blocks in Eq. (Q). For that reason, it is useful to 
make the rotation represented by 



P 



W o 

1 



(33) 



i 

i L I 



T 



i L I 



i 



so that the evolution operator is 



FIG. 5: A plaquette of the AIII model, showing in detail 
the possible scattering processes, and also the effect of the 
rotation specified in Eq. (|j^). Notation is as in Fig.^, and 
the positions of filled circles should be compared in the two 
figures. 



Finally, to put the result in a standard form, we obtain 
H from H using a rotation which takes (a x , a v , er z ) in H 
to (— <7 V , —<7 z ,a x ) in H. We have finally 

H = [{-id x - A x )a x + {-id y - Ay)a y + ma z + V] ® £ z 
+ [B x a x + B y( jy + W] g) E„ . (37) 



o - n CP -i_ ( cos(/3) -sin(/3) 
' " ~ 1 sin(/3) cos(/3) 



U 

o W 



(34) 



In this basis, the bipartite structure of U is clearly also a 
feature of S. The price paid is that the chiral symmetry 
condition, involving Yi x = PE^P -1 , becomes disorder- 
dependent (though not at lowest order in 0) , and for this 
reason we do not use the basis elsewhere in the paper. 
The effect of the rotation can be represented by moving 
some of the points on the network at which wavefunction 
amplitudes are defined, as illustrated in Fig.||. 

In this basis S 2 , like U 2 , is block-diagonal, with sepa- 
rate blocks acting within each of the two sublattices. We 
obtain a Dirac Hamiltonian H by interpreting the part 
of S 2 that acts within one sublattice as exp(— iH), and 
expanding to lowest order in (pi, bi « /3, and m. Using 
the numbering of links around a plaquette specified in 
Fig. ||, we define 



U=(0i+02 + 3 + 04)/2, 

4 -0 2 )/2 and Ay = Ox - 4> 3 )/2 (35) 



as in Ref.|24j. We also define 

lF=(/3 1 +/3 2 +/3 3 +/34)/2, 
B x = ([3 1 -j3 3 )/2 and B y = (fa - fo)/2 . (36) 



As expected, the Hamiltonian H fulfills the chiral sym- 
metry condition, ~E X HY, X = —H. It can be cast in the 
form' of Eq. (U), with 

h = (-id x - A x )a x + {-idy - A y )cr y + mcr z +V. (38) 



Here, A x = A x + iB x , A y 



A. 



iB y and V 



V + iW 



are, respectively, complex random vector and scalar po- 
tentials. If disorder is also introduced in the mass m and 
as an additional term B z a z ® T, y) H is the most gen- 
eral form for a four-component Dirac Hamiltonian with 
randomness in symmetry class AIIIEj. 

Our network model for class BDI can be treated in a 
similar way, but for there to be a continuum limit it is 
necessary to restrict link phases to the value zero, instead 
of allowing either <f>i = or 4>l = We again obtain H 
in the form of Eq. ( |l^ ) , but with 



{-id x -iB x )a x + (-id y 



-iBy)cr v + ma z +iW , (39) 



which is a Dirac Hamiltonian with purely imaginary ran- 
dom vector and scalar potentials and a real mass m. The 
same form is obtained from thp-jcontinuum limit of the 
Hatsugai-Wen-Kohmoto modelj23 which is known to be- 
long to the symmetry class BDI. 
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E. Behaviour for special parameter choices 

The behaviour of the models we have introduced can 
be simplified by making one of several alternative special 
choices for the parameters Pi that characterise scattering 
within the two-component space associated with chiral 
symmetry, or for the parameter a, specifying scattering 
at nodes. We discuss these choices in this subsection, 
restricting ourselves to the model from class AIII, as il- 
lustrated in Fig.|2[ 

1. No chiral scattering 

On setting fti = for every link I, one has two de- 
coupled U(l) network models. With the choices for link 
phases <f>i and node parameters a that are set out in 
Sec.|J, these models are in localised phases for a ^ 7r/4, 
with a localisation length which diverges as a approaches 
tt/4. 



2. Uniform chiral scattering 

There are two alternative and natural ways to pick the 
parameters Pi or 6; without disorder. 

One is to set Pi = Po, a constant. This results in a 
model that is statistically invariant under tt/2 rotations 
of the lattice. Its time evolution operator develops a gap 
in its eigenphase spectrum near E = for /3q 7^ 0. To 
show this, we note from Eq. (|t]) that the eigenphases E a 
of S can be written for constant pi in terms of Pq and 
those of U, which we denote by eg. We find 



cos(£' a ) = cos(/3 ) cos(ea) . 



(40) 



so that there are no eigenphases E a in the interval 
[— Pq,+Po\. Because of the gap, this choice for the Pi 
is of limited interest. 

An alternative choice is to take the hyperbolic angles 
constant, with bi = b in the transfer matrix. This breaks 
the spatial symmetry of the model, because after a tt/2 
rotation of the lattice, the sign of bi is changed in the 
transfer matrix for one half of the links of the system. 
For the transfer matrix before such a rotation, it is clear 
from Eq. ([l4]) that the Lyapunov exponents of the model 
with chiral symmetry, which we denote by v n , are related 
to those of the U(l) model, fi m , by 



= fi m ±b. 



(41) 



We take the fact that the Lyapunov exponents of the two 
models are related in such a simple way as an indication 
that this choice for bi does not lead to behaviour generic 
for the chiral symmetry class. The result is nevertheless 
interesting, because small fluctuations in bi about an av- 
erage value b are likely to restore generic properties of 
the symmetry class without large changes in the values 
of the Lyapunov exponents. In particular, for suitable 



values of b and the node parameter a, the spectrum for 
v n expected from Eq. (|4l]) has a gap around v = 0, indi- 
cating that a localised phase is possible for the network 
model with chiral symmetry. 



3. Strong chiral scattering 

The choice \p{\ = tt/2 for all / breaks the system up 
into a set of disconnected loops, as illustrated in Fig.[|. 
In this limit, possible values for the eigenphases of S are 
E = 0, tx, ±7r/2, , ±a, ±a+7r and ±a±7r/2, where a is the 
scattering angle at nodes. The probability for each eigen- 
value to occur depends on the distribution chosen for the 
signs of Pi . To determine the nature of eigenfunctions as 
this limit is approached, degenerate perturbation theory 
in \P\ — tt/2 is required. At leading order, this generates a 
tight-binding model with nearest neighbour hopping and 
site dilution on a square lattice. 




FIG. 6: The structure of the AIII network model in the lim- 
iting case \/3\ — tt/2. 



4- Decoupled plaquettes 

Setting the node parameter to a = or a = tt/2, the 
system breaks up into decoupled plaquettes. Eigenphases 
for one such plaquette can be calculated analytically, and 
vary continuously with link phases </>;. Eigenstates are 
therefore non-degenerate and localised in this limit. It is 
also interesting to consider behaviour at the edge of the 
system in a localised phase. For a single copy of the U(l) 
network model which is defined on the plane rather than 
a torus, two distinct localised phases arise in the two lim- 
its for a: one with current-carrying edge states and one 
without. Following the construction described in Sec.|J, 
for models with chiral symmetry two oppositely propa- 
gating edge states appear, with scattering between them 
for Pi 7^ 0. Viewed as a strictly one-dimensional system, 
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these edge states are expected to be critical if the aver- 
age (fa) is zero, and localised otherwise. Weak coupling 
between plaquettes, however, is likely to localise the edge 
states for all values of (fii). We therefore do not expect 
an analogue of the quantum Hall plateau transition in 
the models we study, a conclusion that one may also an- 
ticipate from the fact that the models are statistically 
parity-symmetric. 



F. Models with only vector potential disorder 

There has been very extensive interest in Dirac Hamil- 
tonians on two space dimensions, with only vector poten- 
tial disorder. In this subsection we show how lattice ver- 
sions of such problems can be realised, using the network 
model. A similar treatment for tight-binding Hamilixmi- 
ans has been described recently by Motrunich et aZO 

Our starting point is a special choice for the chiral 
scattering angles /?;, derived from a potential function as 
follows. The potential is defined at the points which are 
marked with filled circles on the lattice shown in Fig.||. 
Suppose that, on the network drawn with (say) full lines 
in this figure, the link I runs from point i to point j, 
and let the potential take the values <&i and at these 
points. The hyperbolic angle b is then the gradient of the 
potential, with 



(42) 



Introducing a diagonal matrix <I>, with $/ as entries, 
one finds that the evolution operator hi (6) appearing in 
Eq. (p(]) is a similarity transform of the operator U 



U{b) = e*[/e-* 



Likewise, one has 



U(-b) 



(43) 



(44) 



Also the Hamiltonian defined by Eq.(32) shows a decom- 
position reminiscent of the continuum version 



7i$ = 



e~* 
e* 



where 



Ho — 



Ho 
Ho 



1 1-U 
2H + U 



e" 



,(48) 



(49) 



There are two separate applications of these ideas 
which are of interest. The first is to take the operator 
U to be without disorder, setting all link phases <fri = 0. 
Choosing in addition the node parameter to have its crit- 
ical value a = 7r/4, U has a stationary state ipo which is 
constant in space. Then the wavefunctions Z + and Z_ 
are lattice versions of the eigenstates known for the Dirac 
equation with random vector potential, whose statistical 
properties have been very extensively studied. This is 
supported by the fact that in this case the continuum 
limit of the model Eq. ( |37j ) simplifies to a pair of Dirac 
fcrmions with real vector potential, equivalent to the ones 
described in Ref . [Tt|. 

A second application is to include disorder both via <& 
and in U itself. Taking all link phases to be random, and 
again setting a — 7r/4, ipo is a critical wavefunction for 
the quantum Hall plateau transition. The combinations 



Z± 



ipo will therefore have multifractal fluctuations 



that are determined from an interplay of the distributions 
of ipo an d of <£. It seems very likely that this will result 
in a new type of critical behaviour, but we leave details 
for a future investigation. 



IV. NUMERICAL RESULTS FOR THE 
TWO-DIMENSIONAL AIII MODEL 

A. Introduction 



These facts enables us to determine the stationary 
states for a system with non-zero $ in terms of those for a 
system with $ = 0, if boundary conditions or fine-tuning 
of the disorder realisation allow such states. Specifically, 
suppose that the system without chiral scattering has a 
stationary state so that 



(45) 



Then the system with chiral scattering has stationary 
states that, in the notation of Eq. (pi]), are given by 



Z+ = e^ipo and Z_ = e'^tpo ■ (46) 

A similar simplification applies to the Green function 
at z = 1. The relevant part is the off-diagonal block 
appearing in Eq. (|30|), for which one has 



l + Ue b $1 + U _ $ 

t — e e 

1 - Ue b l-U 



(47) 



In this section we present results from a numerical in- 
vestigation of the two-dimensional AIII network model, 
as defined in Fig. || and Eq. (Q) . Most of our results are for 
systems in which the hyperbolic angles 6/ have a Gaus- 
sian distribution of zero mean and variance g. For these 
systems, we examine behaviour as a function of the node 
parameter a and g, parameterising the latter with 7, de- 
fined by sin(7) = tanh(y / g). Hence we consider the pa- 
rameter space spanned by < a < tt/4 and < 7 < ir/2. 
We also present some data for systems in which the an- 
gles Pi are uniformly distributed between and 2n. To 
set the scene, we first summarise what can be anticipated 
for behaviour of the model, from a consideration of limit- 
ing values for the parameters a and 7 and from existing 
analytical results. 

We start by recalling that at 7 = the system con- 
sists of two decoupled copies of the U(l) network model. 
States in each copy are localised for all a in the range 
< a < 7r/4, with a localisation length that diverges at 
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/A follow from 
the renormal- 



(50) 



a = 7r/4. One expects the localised phase to be stable 
against small perturbations, so that it will extend over 
a region in parameter space with < a < tt/A and 7 
small. Conversely, critical states in the two U(l) models 
at a — tt/A are easily mixed by non-zero chiral coupling 
bi , so that it is reasonable to guess-that the critical phase 
predicted by Gade and Wegnertffl will be found in the 
region with a close to tt/A and 7 > 

Properties along the (critical) line 
studies of the non-linear sigma model 
isation group equations 

d\A _ ^_ 
~dT ~~ 16 

are are uedEl to be exact to all order of perturbation the- 
ory. Note that when a is close to but not exactly equal 
to 7r/4, a non-zero mean mass for the underlying ran- 
dom Dirac fermions is generated. In that case, Eq. j50| ) 
should be extended to incorporate it, which is beyond 
the scope of the present work. This non-zero mass leads 
to localisation when large enough on a scale set by 7. 

The coupling constant A parameterises the conductiv- 
ity er, with (7 = l/(87rA 2 ). Since A does not flow under 
renormalisation, the critical phase is described by a line 
of fixed points and may have a range of values for its 
conductivity. Flow of the second coupling constant, Xa, 
determines properties at non-zero energy, giving in par- 
ticular for the jdipsity of states p(E) at small energy E 
the scaling lawBu 

En 



p{E) ~ exp (- Kv /ln(£o/£)) , 



(51) 



where Eq is a microscopic energy scale. Using the results 
of Ref.|l|, k = 4\/2/A 2 . As emphasised in Ref.[l7| such 
behaviour for p(E) sets in only below a crossover energy 
scale E c : from Ref. [l6| we estimate this scale to be given 
by 



E c = E cxp (-acr 2 ) , 



(52) 



where a = (32^/2n) 2 . We note also that Motrunich et 
al have argued that the power law in the exponent of 
Eq. (|5l]) is modified at very small energy scalesjill and 
that their results liaye been reproduced within a field- 
theoretic approach.^ 

Since the derivations of these results depend on a num- 
ber of technical assumptions, some checks are very desir- 
able. Unfortunately, computational work is challenging 
because it is hard to reach scales smaller than E c , and 
at intermediate energies the same analysis predicts an 
approximately power-law dependence of the density of 
states on energy, with an effective exponent which de- 
pends on disorder strength. 

In our numerical work testing these expectations, we 
use calculations of two types, both applied to quasi-one 
dimensional systems of width M plaquettes and length L 
plaquettes (so that the time evolution operator S appear- 
ing in Eq. (M) is a AML x AML matrix, and the trans- 
fer matrix T appearing in Eq. ( [l4] ) is a AM x AM ma- 
trix). For such samples, we find the localisation length 



£m as the inverse of the smallest positive Lyapunov ex- 
ponent for the transfer matrix, calculated using the stan- 
dard approach.^ Applied directly to the transfer matrix 
as formulated in Sec.|l| the resulting localisation length 
is that for eigenstates of the time evolution operator at 
eigenphase zero; by associating an extra phase e lE with 
every link (breaking chiral symmetry) , we are also able to 
find the localisation length for states with non-zero eigen- 
phase. Separately, we determine the density of states for 
eigenphascs of the time evolution, operator S, applying a 
version of the recursion method,E3 adapted to treat a uni- 
tary operator in place of the Hermitian one for which the 
approach was originally formulated. We use the terms 
eigenphase and energy interchangeably in describing our 
results. 



B. Phase diagram 

We determine a phase diagram for the model in the 
a, 7 plane from a study of the zero-energy localisation 
length in quasi-one dimensional geometry, and its de- 
pendence on M, (We do not apply this approach to 
models with uniformly distributed since our trans- 
fer matrix calculations develop numerical instabilities if 
|cos(/3;)| <C 1.) The quantity of central interest is the 
ratio £m/M. In a localised phase, £m approaches a finite 
limit £, the bulk localisation length, for large M, and so 
the ratio £,m/M decreases with increasing M, varying as 
£/M for sufficiently large M. By contrast, in a critical 
phase, the ratio approaches for large M a finite constant 
whose value, can be identified with the conductivity of 
the modelO 

Representative data for £,m/M at M = 8, 16,32 and 
64 are shown as a function of 7, in Fig.^ for a = 7n/32, 
and in Fig.H for a — ir/8. In both cases, two regimes 
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FIG. 7: £m/M at a = 7ir/32 as a function of 7. 

of behaviour are observed: for small 7, £,m/M decreases 
with increasing M, indicating a localised phase, while 
for 7 close to 7r/2, £,m/M is independent of M for the 
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FIG. 8: £m/M at a = 7r/8 as a function of 7. 



system widths studied, suggesting a critical phase. By 
identifying the value of 7 which divides the two regimes 
and studying its dependence on a, we arrive at the phase 
diagram displayed in Fig|| 
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FIG. 9: Phase diagram of the AIII network model in the (a, 7) 
plane. 

Properties in the critical phase are illustrated in more 
detail by the behaviour of the ratio £m/M as a function 
of 7 with a = 7r/4, shown in Fig.^j. On this line in pa- 
rameter space, behaviour at all points is critical rather 
than localised, and the data show a finite limiting value 
a for the ratio at large M, with a dependent on 7: such 
behaviour is expected from Eq. ( |50| ) and the identifica- 
tion of the coupling constant A as being a function only 
of conductivity, which is given in turn by £m/M. This 
figure also illustrates crossover of behaviour with increas- 
ing M at small 7, from that of the U(l) model to that of 
the AIII model. The crossover is responsible for discon- 



tinuous variation of a with 7 at 7 = 0. 
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FIG. 10: £m/M at a — tt/4 as a function of the 7. Inset: the 
limiting value, a, of this ratio for large M 



We show in Fig. [ll] the energy dependence of the locali- 
sation length at a point in the localised phase (left panel) 
and at a point in the critical phase (right panel), with 
data at L = 2.5 x 10 5 for a range of system widths M, as 
well as an extrapolation to the two-dimensional limit. It 
is clear from the data that the localisation length in the 
critical phase decreases very rapidly with E away from 
E = 0. 




icf 9 icr 7 lcr 5 ltr 3 kt 9 icr 7 icr 5 icr 3 



FIG. 11: Localisation length £m at a = it/8, g — 1 equivalent 
to 7 = 0.55-7r/2 (left) and at a — n/4, g = 9 equivalent to 
7 = 0.936-7T/2 (right) as a function of E in a log-log scale 
for increasing transverse size M=2,4,8,16,24 and extrapolated 
to Af = co, indicated by the symbols (), A,0,\j,Q and * 
respectively. 



12 



C. Density of states 

Our interest in the (average) eigenphase density of the 
time evolution operator is focussed mainly on behaviour 
near zero energy, where we expect in the critical phase 
singular behaviour approaching that of Eq. (|5l]), and in 
the localised phase either a density vanishing quadrati- 
cally with energy, asJbr the random matrix ensemble in 
this symmetry class ,□ or singularities arising from Grif- 
fiths strings, as discussed in Ref. 



the large values of z reached here indicate close approach 
to the limiting low-energy behaviour. 
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Before examining 
the small energy region in detail, we give an overview of 
behaviour for the entire range of eigenphases, shown in 
Fig. |l2] for a = ir/4 and uniformly distributed fti. (Simi- 
lar results are obtaine d with Gaussian bi at g ~ 1.) From 
the discussion of Sec. [II A| , we need to display data only 
over the range < E < ir/2. These results were obtained 
for a system size of M = 16 and L = 10 5 , which is large 
enough for p(E) to be self-averaging at the scale used 
here. 
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FIG. 12: Density p(E) of eigenphases E for a = 7r/4 and /3; 
uniformly distributed. 

A divergence of p(E) as E — > is apparent in Fig.jl^ , 
which we now investigate in more detail. In Fig.O we 
show p(E) as a function of E, with logarithmic scales for 
both axes, for a sequence of models, with a — 7r/4 and 
bi Gaussian distributed in each case. Members of the 
sequence are chosen to have successively smaller values of 
a, and hence larger values for the energy scale E c , since 
by increasing 7 we decrease cr, as is evident from the 
inset to Fig. ^. The data were obtained using systems of 
width M = 16, and lengths in the range L = 10 5 to L = 
10 9 , in order to achieve satisfactory self-averaging. Over 
the energy range accessible, the variation of p{E) with 
E is approximately power-law, and can be characterised 
following Ref. [l?] using an effective dynamical exponent 
z, with 

p(E)~ \E\- 1+2 ' Z . (53) 

As one approaches the asymptotic low-energy behaviour, 
the effective exponent z increases, ultimately to infinity: 



FIG. 13: Density of states p(E) as a function of E for small 
E, with logarithmic scales on both axes. Three cases of in- 
creasing disorder are plotted (with symbols A for g = 400, 
for g — 25 and V f° r 9 = 9)- The related dynamical expo- 
nents z in this range of energy are respectively z — 26.0 ±2.0, 
z — 8.0 ± 0.1 and 5.3 ± 0.1. Symbols □ show the case f3i 
uniformly distributed in [0, 2-k] with z = 3.02 ± 0.03. 

By repeating such calculations for a range of M values, 
we have checked that these results are not influenced by 
finite system width system. This is shown in Fig. [l4|, in 
which the behaviour with uniformly distributed 0i is also 
illustrated. The rather rapid convergence to the large M 
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FIG. 14: Density p(E) of eigenphases E for a = 7r/4, g — 100 
and for increasing transverse sizes M=4,8,16 and 24. The 
inset corresponds to the same set of system sizes but with /3i 
uniformly distributed in [0, 2ir], 

limit is a consequence of the small values of the locali- 
sation length for energies away from E = 0, illustrated 
above. 
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Finally, we examine the evolution of the density of 
states as one moves from the critical phase into the lo- 
calised phase. Results are presented in Fig. [ll| At the 
broad scale to which these calculations are restricted, be- 
haviour in the localised phase is consistent with the power 
law of Eq. (|5^), with a power —1 + 2/z which, deep in 
the localised phase, is positive and increases as the local- 
isation length decreases. Such behaviour was proposed 
as generic for localised systems with chiral symmetry in 
Ref.0 
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V. SUMMARY 



In summary, we have introduced network models which 
are realisations of the chiral symmetry classes. We have 
argued that they are interesting from several points of 
view. They are useful as a starting point for numerical 
studies of these symmetry classes, offering access in the 
two-dimensional case we have examined to both critical 
and localised phases, and displaying band-centre singu- 
larities in the critical phase which approach quite closely 
the expected asymptotic form. The models also have 
a striking connection to network models without chiral 
symmetry, but with absorption and amplification. More- 
over, by imposing constraints on the disorder, they serve 
as lattice versions of problems with randomness enter- 
ing only through a vector potential. It seems likely that 
disorder of this kind can generate types of critical be- 
haviour different from those known previously for locali- 
sation problems in two dimensions. 
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